#!/usr/bin/perl -w
use strict;

# an example of what we will parse is
# >ARB_00831 |  ARTHRODERMA BENHAMIAE CBS 112371 NULL (410 AA) [Cleavage at C-17 of 50% quality]

# This script will parse GPI-SOM output and find the cleaved protein length
# The output will just print the size of protein and the size of the protein
# after it is cleaved

use Bio::DB::Fasta; # database to pull in the full sequence

# script to take 2 arguments
# 1. Protein database
# 2. GPI-SOM output
my ($db,$gpioutput) = @ARGV;

# run it like: perl parse_GPISOM.pl proteindbfile gpisomfile

my $dbh = Bio::DB::Fasta->new($db); # index the FASTA database

open(my $fh => $gpioutput) || die "Cannot open $gpioutput: $!\n";

print join("\t", qw(PROTEIN_NAME LENGTH C_POSITION 
		    CLEAVED_LENGTH PROBABILITY)), "\n";
# process each line in the file
while(<$fh>) {
    if( /^>(\S+)\s+.+\[Cleavage\s+at\s+C\-(\d+)\s+of\s+(\d+)\%/ ) {
	my ($protein,$c_position, $prob) = ($1,$2,$3);
	my $seq = $dbh->get_Seq_by_acc($protein);
        if( ! defined $seq ) {
	   warn("cannot find sequence $protein\n");
	   next;
	}
	my $cleaved_seq = $seq->trunc(1, 
				      $seq->length - $c_position);
	print STDERR "ALL    ", $seq->seq, "\n";
	print STDERR "CLEAVE ", $cleaved_seq->seq, "\n";
	print join("\t", 
		   $protein,
		   $seq->length,
		   $c_position,
		   $cleaved_seq->length,
		   $prob,
		   ),"\n";
    } else { 
#	warn "ignoring line that doesn't match: $_";
    }
}


